home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Windows Expert
/
Windows Expert.iso
/
desktop
/
profft.zip
/
SOURCE
/
FFT.CPP
< prev
next >
Wrap
C/C++ Source or Header
|
1992-04-27
|
4KB
|
150 lines
/****************************************************************************
FFT.CPP Denne filen inneholder (i denne rekkef°lgen):
*****************************************************************************
#define M_PI 3.14159265358979323846
int sincos_length = 0;
int newpow = 0;
prcomplex *sincos_tab;
int sincos(int length)
int fft(prcomplex *z1, prcomplex *z2, int n)
****************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#pragma hdrstop
#ifndef __PRCOMP_H
#include "prcomp.h"
#endif
// Definerer M_PI som brukes ved oppretting av sinus cosinus tabellen
#define M_PI 3.14159265358979323846
// Husker lengden pσ tabellen
int sincos_length = 0;
int newpow = 0;
// Peker til tabellen
prcomplex *sincos_tab;
int sincos(int length)
/****************************************************************************
Denne lager en tabell for raskt oppslag pσ komplekse sinus og cosinusverdier
som brukes av fft algortimen. Denne er sakset fra et dokument som ble delt
ut i forbindelse med oppdraget og stammer fra Universitetet i Oslo.
int length Lengden pσ tabellen.
Returnerer: En toerpotens av st°rrelsen (?!?).
Kodet av: MK 02.04.92 Modifisert for vσre komplekse talltyper.
*****************************************************************************/
{
int i, totlen;
// Hvis tabellen er opprettet allerede
if (sincos_length == length) return(newpow);
// Hvis lengden er ulik 0 men ikke lik lengden sσ fjern tabellen
// fra minnet.
if (sincos_length != 0) free(sincos_tab);
// Sett riktige startverdier
sincos_length = 0;
i = 1; newpow = 0;
while (i < length)
{
i <<= 1;
newpow++;
}
if (i != length)
newpow = 0;
totlen = length << 1;
// Alloker minne til tabellen
sincos_tab = (prcomplex*) malloc(sizeof(prcomplex)*(totlen+1));
if (sincos_tab==NULL)
return(-1);
// Ta vare pσ lengden av tabellen til senere
sincos_length = length;
for (i=0; i<=totlen; i++)
{
// Fyll inn de kalkulerte verdiene i tabellen.
sincos_tab[i].re = (real)cos((double) i / (double) length * M_PI);
sincos_tab[i].im = (real)-sin((double) i / (double) length * M_PI);
}
return(newpow);
}
int fft(prcomplex *z1, prcomplex *z2, int n)
/****************************************************************************
Dette er en rask FFT algoritme sakset fra et dokument som ble delt
ut i forbindelse med oppdraget og stammer fra Universitetet i Oslo. Den
er derfor ikke dokumentert i noen sµrlig grad.
prcomplex *z1 Peker til f°rste arbeidsbuffer. Mσ vµre ferdig allokert
og utfyllt ved kall.
prcomplex *z2 Peker til andre arbeidsbuffer. Mσ vµre ferdig allokert.
Returnerer: 1 Hvis f°rste arbeidsbuffer inneholder transformen.
2 Hvis andre arbeidsbuffer inneholder transformen.
0 Hvis transformen feiler.
Kodet av: MK 02.04.92 Modifisert for vσre komplekse talltyper.
MK 03.04.92 Fjernet un°dvendig kode.
*****************************************************************************/
{
int power, inda, nh, dir;
int zstep, inzee, ind, ia;
prcomplex *ar1, *ar2;
prcomplex arg2;
prcomplex *zind, *fpt, *tpt, *pt;
real fn;
// Hvis tabellen ikke kunne lages.
if (sincos(n)<=0) return(0);
nh = n / 2;
power = dir = 1;
ar1 = z1;
ar2 = z2;
for (zstep=n;zstep!=1;zstep>>=1)
{
zind = &sincos_tab[n];
ia = power;
while (ia--)
{
zind = &zind[-zstep];
fpt = &ar1[inda=nh+ia-power];
tpt = &ar2[2*inda-ia];
while (inda >= 0)
{
pt = &fpt[nh];
arg2.re = zind->re * pt->re - zind->im * pt->im;
arg2.im = zind->re * pt->im + zind->im * pt->re;
pt = &tpt[power];
tpt->re = fpt->re + arg2.re;
tpt->im = fpt->im + arg2.im;
pt->re = fpt->re - arg2.re;
pt->im = fpt->im - arg2.im;
fpt = &fpt[-power];
tpt = &tpt[-power];
tpt = &tpt[-power];
inda -= power;
}
}
pt = ar1;
ar1 = ar2;
ar2 = pt;
dir = 3 - dir;
power <<= 1;
}
return(dir);
}